Optimal non-pharmaceutical pandemic response strategies depend critically on time horizons and costs

The COVID-19 pandemic has called for swift action from local governments, which have instated non-pharmaceutical interventions (NPIs) to curb the spread of the disease. The swift implementation of social distancing policies has raised questions about the costs and benefits of strategies that either aim to keep cases as low as possible (suppression) or aim to reach herd immunity quickly (mitigation) to tackle the COVID-19 pandemic. While curbing COVID-19 required blunt instruments, it is unclear whether a less-transmissible and less-deadly emerging pathogen would justify the same response. This paper illuminates this question using a parsimonious transmission model by formulating the social distancing lives vs. livelihoods dilemma as a boundary value problem using calculus of variations. In this setup, society balances the costs and benefits of social distancing contingent on the costs of reducing transmission relative to the burden imposed by the disease. We consider both single-objective and multi-objective formulations of the problem. To the best of our knowledge, our approach is distinct in the sense that strategies emerge from the problem structure rather than being imposed a priori. We find that the relative time-horizon of the pandemic (i.e., the time it takes to develop effective vaccines and treatments) and the relative cost of social distancing influence the choice of the optimal policy. Unsurprisingly, we find that the appropriate policy response depends on these two factors. We discuss the conditions under which each policy archetype (suppression vs. mitigation) appears to be the most appropriate.


The SIR Model
We use a standard Susceptible, Infectious, Recovered (SIR) model of disease transmission. The standard SIR model describes disease transmission in a well-mixed system in which an interaction between any two individuals is equally likely. Our approach implements the non-dimensional form of the standard SIR model because we couple it to the dimensionless forms of cost functions to find solutions to the social distancing problems. We start by providing a derivation of the nondimensional SIR form starting from its standard form.
The differential equations describing the dynamics of the susceptible, infectious, and recovered individuals in the population are: where S is the number of the susceptible, I is the number of the infectious, R is the number of the recovered, and N is the total number of individuals in the population. The parameters β and γ respectively represent the transmission and recovery rate and have units of inverse time.

Non-Dimensionalization
As is typically done, we define the dimensionless unit of time τ , which is: In other words, in the dimensionless SIR model, τ = 1 represents the average duration an individual is infectious for. We further recast the state variables in their dimensionless form where The dimensionless differential equations for the population dynamics become: By defining the reproduction number R 0 = β γ , these equations can be formulated as

Final Epidemic Size
Without mitigation measures, and under the assumption that the initial conditions for the proportion of the population infected is very small (i.e., i(τ = 0) ≪ 1), and nearly the entire population is initially susceptible (i.e., s(τ = 0) ≈ 1), the proportion of the population that remains susceptible as τ → ∞ is given by the solution to the equation [Diekmann(2013)] and the final epidemic size is 1 − s ∞ .

The Effective Reproduction Number and Herd Immunity
The effective reproduction number R τ describes the average number of infections produced by each infected individual in the population at time τ . In the standard SIR model without mitigation measures, R τ = sR 0 .
Herd immunity occurs when the number of those susceptible in the population is reduced either by natural infection or vaccination such that the effective reproduction number is less than 1. Therefore, the maximum proportion of susceptible individuals in a population with herd immunity is

Time-Dependent Transmission
In classic SIR models, β and R 0 are constant parameters. This paper considers the case where the transmission rate, β, is dynamic and influenced by behavioral changes. We will use the notation that β is a dynamic variable and β 0 is the transmission rate without any behavioral mitigation. Similarly, we define R D (τ ) as the dynamic reproduction number, which is a function of social distancing behavior. R 0 is the reproduction number in the absence of any behavioral mitigation. In a dynamic framework, the dimensionless SIR model, Equations 5 become: In the framework with a dynamic transmission rate, the effective reproduction number becomes:

Cost Functions
We want to find the social distancing policy that minimizes a total cost, which includes the cost of social distancing and the cost of infections over a pre-specified time horizon ranging from time t = 0 to t = t final . The social distancing policy is described by β(t), the transmission rate, that depends on social distancing behavior at time t. We further define β 0 to be the infection rate for the disease in the absence of any behavioral change. Therefore, β(t) is restricted to the interval (0, β 0 ] The total cost is: where H sd (t) is the cost of social distancing per unit time, and H infect is the cost of infections. As the transmission rate β(t) decreases, the cost of social distancing per unit time, H sd (t), increases, but the cost of infections per unit time H infect (t) decreases.

Cost of Infections
The number of new infections per unit time is βSI N , which comes from the SIR equations. If we define D to be the average cost of infection per infected individual. Then, the cost per unit time of infections is: Note that the total cost of infections is:

Cost of social distancing
We assume that the total cost of social distancing is proportional to the size of the population N , and a cost parameter C, which has units cost per person per unit time. Therefore, we define the cost per unit time of social distancing to be: Here, g(x) is the function that defines the relationship between the relative cost of social distancing and relative reduction in the transmission parameter. The function g(x) should have the following properties: 1. g(x) is a monotonically decreasing function on the interval [0, 1]. The theoretical setup of the problem determines this property. The cost of social distancing should increase as the transmission rate is further decreased.
2. lim x→0 + g(x) = ∞. In theoretical terms, it is not possible to stop all transmission completely; therefore, the cost of decreasing transmission rates to zero should be infinite. From a practical perspective, this restriction will prevent optimal solutions to β(t) from passing through β(t) = 0, and if the initial condition for β(t) is positive, the solution will remain positive.
3. d dx g(x)| x=1 = 0. Theoretically, if there is no cost of infections, the optimal β(t) = β 0 . Therefore, the cost of social distancing should be minimized when β(t) β 0 = 1. From a practical perspective, combined with condition 2, this condition ensures that solutions to β(t) will be bounded to the interval (0, β 0 ] if correctly initialized to the interval. 4. g(1)=0. When β(t) β 0 = 1, there are no behavioral changes so the cost of social distancing should be zero.
Based on these three required properties, we choose the function to have the form and thus, substituting this into equation 15, the cost per unit time of social distancing is

Dimensionless Cost Equations
We now derive the dimensionless form of the cost equation. We define a dimensionless cost of social distancing, c to be: Effectively, c compares the cost of social distancing to the cost of infection. Recall that C has dimensions of cost per individual per unit time, D has dimensions of cost per individual, and γ has units of inverse time. Therefore, c is dimensionless. When c is small, people will more willingly social-distance as they feel that the cost of becoming ill outweighs the benefits of mixing socially. Their preference favors their lives and health compared to livelihoods and non-health-related wellbeing. When c is large, people's preferences are switched. We define the dimensionless cost: We can make the following substitutions in all equations: Then, the dimensionless cost of infection, Equation 13, becomes: The dimensionless cost of social distancing, Equation 15, becomes: We denote the cost function as h(s, i, R D ) which can be expressed as: By setting τ final = γT final , the total cost in its dimensionless form is
The functions y j [Weinstock (1974)] that satisfy the following conditions are extrema of the integral in Equation 25: where y ′ j indicated the derivative dy j dτ . Our main objective is to solve for R D that minimizes Equation 23. Because the integrand is a function not only of R D , but of s and i as well, we introduce the functional Lagrange multipliers λ i and λ s to ensure that the conditions of the differential equations governing the infection dynamics (Equations 5) are met. Then, based on the constraints provided by the SIR dynamics, the full Lagrangian becomes: Note that we do not need to include a differential equation constraint for r since the equations for s and i ensure that r = 1 − s − i.
We can now apply Equation 26 to all y j ∈ {s, i, R D , λ i , λ s }. y j = s yields: y j = i yields: y j = R D yields: Note that applying y j = λ s and y j = λ i to Equation 26 returns the constraint differential equations for s and i from Equation 5. Note that applying y j = λ s and y j = λ i to Equation 26 returns the constraint differential equations for s and i from Equation 5. The differential equations can then be solved as a two-point boundary value problem. For all cases presented in this paper, we set i(0) = i 0 , s(0) = 1 − i 0 , λ s (τ final ) = 0, and λ i (τ final ) = 0.
[rv]the τ is not being displayed correctly -anyone know why? We set λ s and λ i to be equal to zero at the end point of the interval because s(τ final ) and i(τ final ) are unconstrained. We solved the 2-point boundary value problems in R using the bvpSolve package. [Mazzia et al.(2014)Mazzia, Cash, and Soetaert]

Multi-Objective Formulation
So far, we have framed our problem as a single-objective optimization. In other words, we optimized the sum of a social distancing and infection costs, where the parameter c defines how a level of social distancing should be weighed against a level of infections. In a Pareto-optimization formulation, we find the solutions R D (τ ) for which the total social distancing cannot be decreased without increasing infections and infections cannot be decreased without increasing the social distancing. Mathematically, the total social distancing, which we call k SD is: and the cost of infections, k infect is given by: We recover the total cost defined in Eq. 24 though the equation We know that k infect is the total number of new infections that occur between time τ = 0 and time τ = τ final , which means that k infect can be defined as: The only term in Eq. 34 that depends on R D (τ ) is −s(τ final ). We can therefore formulate the set of Pareto-optimal solutions as those that minimize k SD for each possible value of s(τ final ).
The Lagrangian in this formulation is then given by: The variables µ s and µ i are Lagrange multipliers in the multi-objective formulation.
We can now apply equation 26 to the variables in Equation 35. Applying y j = s yields: y j = i yields: y j = R D yields: Now, we show that the formulations are equivalent if λ s = cµ s − 1 and λ i = cµ i Starting from Equation 28, we have: Which becomes: Starting from Equation 29, we have: which becomes Starting from Equation 30, we have: which becomes In the alternative formulation, s(τ final ) = s final . In this formulation, because s is fixed at the final endpoint, only µ i = 0. We know that We know that λ s (τ final ) = 0. We can therefore calculate c if we solve the equations in our second formulation as follows: µ s (τ final ) = 1 c (48)

Contributions
• Sarah Nowak (Ph.D.) led the effort in conceptualizing and formulating the model structure and led the overall modeling effort, including model implementation. Sarah is the corresponding author of this study and can be contacted at: sarah.nowak@med.uvm.edu • Pedro Nascimento de Lima (MSc) contributed to the model analysis and interpretation and drafting and revising the manuscript.
• Raffaele Vardavas (Ph.D.) contributed to the model formulation, implementation, and interpretation of model outputs and drafting and revising the manuscript.